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16.1  Abstract 


In  the  study  of  mixing  and  reaction  in  turbulent  flows,  there  are  several  phenomena  that 
can  be  usefully  described  in  terms  of  surfaces.  Examples  are  turbulent  flames  and  the 
turbulent  mixing  of  different  liquids.  The  most  fundamental  type  of  surface  is  the  mate¬ 
rial  surface  which,  by  definition,  moves  with  the  fluid.  Because  of  the  fluid's  turbulent 
motion  and  deformation,  the  surface  is  continually  stretched  and  bent.  In  this  study 
numerical  simulations  have  been  performed  to  understand  and  to  quantify  these  proc¬ 
esses. 

A  pseudo -spectral  method  is  used  to  solve  the  Navier-Stokes  equations  which  govern 
the  motion  of  the  fluid.  These  equations  are  solved  on  a  1283  grid  for  the  simplest  pos¬ 
sible  turbulent  flow  —  statistically  stationary,  homogeneous,  isotropic  turbulence.  As  the 
results  show,  the  direct  numerical  representation  of  a  material  surface  is  not  feasible:  for 
the  surface  area  grows  exponentially  (by  a  factor  oflO17  over  the  duration  of  the  simu¬ 
lations);  and  radii  of  curvature  less  than  a  millionth  of  the  grid  spacing  arise.  Instead  an 
indirect  method  is  used  in  which  ensembles  (4-8,000)  of  infinitesimal  surface  elements  are 
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followed.  Statistics  of  interest  are  obtained  from  the  stretching  and  curvatures  of  these 
elements. 

For  the  first  time,  the  mean  rate  of  stretching  has  been  determined.  It  is  found  that  the 
surface  area  doubles  every  2  1/2  Kolmogorov  time  scales.  (The  Kolmogorov  time  scale 
is  the  smallest  physical  time  scale  in  turbulence.)  While  this  is  certainly  rapid  growth,  it  is 
only  40%  of  theoretical  estimates,  for  reasons  that  are  explained.  Hitherto,  little  has 
been  known  about  the  curvature  of  material  surfaces.  The  results  show  that  extremely 
small  radii  of  curvatures  arise,  as  small  as  10-8  of  a  Kolmogorov  length  scale  (the 
smallest  turbulent  scale).  These  highly  curved  elements  are  found  to  be  almost  perfectly 
cylindrical  in  shape.  Many  other  more  refined  statistics  have  been  obtained. 

The  numerical  simulations  were  performed  on  an  IBM  3090 -600 E,  with  full  exploitation 
of  its  vector,  parallel  and  large-memory  facilities.  A  typical  run  requires  a  total  of  80 
CPU  hours,  but  can  be  completed  in  20  hours  because  all  six  processors  are  used  in 
parallel. 
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Stretching  and  Bending  of  Material  Surfaces  in 
Turbulence 

16.2  Introduction 


In  the  work  described  here,  supercomputer  calculations  are  used  to  answer  funda¬ 
mental  questions  of  long  standing  about  the  effects  of  turbulence  on  material  sur¬ 
faces. 

Turbulence  —  because  of  its  importance  in  the  atmosphere,  the  oceans,  engi¬ 
neering  equipment  and  elsewhere  —  has  been  the  subject  of  theoretical  study  for 
over  75  years.  In  many  areas,  including  the  study  of  material  surfaces,  the  theory 
is  built  on  hypotheses  that  have  not  been  tested  because  of  insuperable  exper¬ 
imental  difficulties.  Now,  however,  by  simulating  turbulence  using  a  supercom¬ 
puter  it  is  possible  to  test  these  hypotheses.  But  much  more:  the  simulations 
provide  fresh  insights  into  the  basic  physical  processes  of  turbulence. 

A  material  surface  is  defined  by  its  initial  condition  (e.g.  a  specified  plane  at 
t  =  0),  and  by  the  condition  that  every  point  on  the  surface  moves  with  the  local 
fluid  velocity.  We  now  provide  two  examples  to  illustrate  the  physical  signif¬ 
icance  of  these  surfaces. 

Figure  147  shows  a  sketch  of  the  mixing  of  two  bodies  of  water  (A  and  B)  in  a 
closed  vessel.  Initially  (t  =  0)  A  contains  a  trace  solute  of  concentration  cj>  =  <£0, 
while  B  is  pure  water  (4>  =  0).  We  consider  the  material  surface  that  is  initially 
coincident  with  the  interface  between  A  and  B.  The  water  is  set  in  turbulent 
motion  and  as  a  result  the  material  surface  is  convected,  stretched  and  bent.  In 
the  first  stages  of  mixing,  the  concentration  is  uniform  everywhere  except  in  the 
immediate  vicinity  of  the  material  surface  where  there  is  a  thin  diffusive  layer. 
Locally,  the  behavior  of  the  diffusive  layer  is  (to  an  excellent  approximation)  the 
same  as  a  plane  layer  in  a  uniform  strain  field.  Consequently,  the  mixing  process 
can  be  completely  analyzed  in  terms  of  the  statistics  of  straining  on  the  material 
surface  [1,  2,  3].  (At  later  times  the  material  surface  folds  over,  and  the  analysis 
breaks  down  once  the  distance  between  folds  is  comparable  to  the  diffusive-layer 
thickness.) 
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(a)  (b) 

Figure  147.  Mixing  Between  Two  Bodies  of  Water,  sketch  of  the  mixing  between  two 
bodies  of  water. 

a)  The  material  surface  that  initially  (t  =  0)  separates  A  from  B. 

b) Normalized  concentration  profiles  normal  to  the  surface.  (Note  6  <  <  L.) 

The  second  example  is  an  idealization  of  a  turbulent  premixed  flame  (such  as 
that  in  a  spark-ignition  automobile  engine).  Under  appropriate  conditions  [4], 
there  is  a  thin  flame  sheet  that  forms  a  connected  but  highly  wrinkled  surface  that 
separates  the  reactants  from  the  products  (see  Figure  148).  This  flame  surface  is 
convected,  bent  and  strained  by  the  turbulence,  and  propagates  (relative  to  the 
fluid)  at  a  speed  w.  If  the  propagation  speed  w  is  small  compared  to  the  turbu¬ 
lent  velocity  scales  (in  a  way  made  precise  in  [3],  then  the  flame  surface  behaves 
like  a  material  surface.  For  this  case  the  statistics  of  interest  are,  again,  the 
straining  on  the  material  surface  and  also  its  curvature  —  for  this  affects  the  prop¬ 
agation  speed  w  [5,  6]. 

The  straining  on  material  surfaces  was  first  and  most  comprehensively  studied 
by  Batchelor  [1,7]  over  30  years  ago.  He  introduced  two  conjectures: 
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Reactants 

Figure  148.  Development  of  a  Turbulent  Flame  Sheet,  sketch  of  a  turbulent  flame 
sheet  developing  from  an  initially  spherical  kernel. 


1.  straining  is  persistent  (i.e.  the  time  scale  of  change  of  the  strain  rate  is  large 
compared  to  the  time  scale  of  strain  itself); 

and 

2.  the  material  surface  becomes  aligned  with  the  principal  axes  corresponding  to 
the  two  greatest  principal  strain  rates. 

Subsequent  workers  (see  Monin  and  Yaglom  [8]  for  a  review)  have  accepted 
these  conjectures  and  Corrsin  [9]  provides  additional  quantitative  arguments  for 
the  persistence  of  strain. 

The  curvature  of  material  surfaces  has  not  been  examined  in  as  much  detail,  but 
it  is  often  assumed  (e.g.  Klimov  [10])  that 

3.  The  radii  of  curvature  of  a  material  surface  are  no  smaller  than  the  smallest 
turbulence  scales. 


Because  of  insuperable  experimental  difficulties,  these  conjectures  could  not  be 
tested  experimentally.  Our  results  show  that  none  of  them  is  correct. 

The  objective  of  the  work  described  is  to  use  simulations  of  turbulence  to  study 
the  processes  affecting  material  surfaces,  in  particular  to  characterize  the  statistics 
of  surface  straining  and  curvature. 
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The  computational  simulation  of  the  turbulence  is  described  in  the  next  section, 
while  the  algorithm  used  to  determine  surface  properties  is  described  in  16.4, 
“Evolution  of  Material  Surfaces.”  The  results  are  given  in  16.5,  “Results  and 
Discussion,”  and  the  paper  closes  with  conclusions  in  16.6,  “Conclusions.” 

Space  does  not  permit  all  the  details  to  be  presented  here:  more  information 
can  be  found  in  previous  works  of  the  authors  [2,  3,  11,  12,  13]. 


16.3  Turbulence  Simulations 

16.3.1  Pseudo-Spectral  Method 


The  governing  equations  of  fluid  motion  —  the  Navier-Stokes  equations  —  are 
solved  for  the  simplest  possible  turbulent  flow:  homogeneous,  isotropic, 
statistically-stationary  turbulence.  Since,  according  to  the  Kolmogorov  1941 
hypotheses  [8],  the  small  scales  of  turbulence  are  universal,  the  study  of  this  sim¬ 
plest  of  turbulent  flows  has  broad  significance. 

The  Navier-Stokes  equations  are  solved  numerically  using  the  pseudo-spectral 
method  developed  by  Rogallo  [14].  The  time-dependent  Eulerian  velocity  field, 
u(x,  t),  is  represented  on  an  equispaced  grid  of  N3  grid  points  which  form  a  cubic 
computational  domain  of  length  Lo.  In  the  simulations  reported  N  =  64  or  N 
=  128.  The  velocity  field  is  continued  periodically  (i.e.,  u(x  +  zLo,  t)  =  u(x,  t), 
where  z  is  any  integer  vector),  and  consequently  u(x,  t)  has  a  finite  Fourier  repre¬ 
sentation.  There  are  N3  corresponding  discrete  nodes  in  wavenumber  space.  Let 
k  be  the  wavenumber  vector  at  a  given  node,  and  k  be  its  magnitude.  The  lowest 
nonzero  wavenumber,  denoted  by  k0,  is  2^/Lo.  The  components  of  k  are  integer 
multiples  of  k0,  ranging  from  (l-N/2)k0  to  (N/2)k0.  We  denote  by  u(k,  t)  the 
complex  Fourier  (wavenumber)  coefficients  of  the  velocity  at  time  t:  i.e.,  u(k,  t) 
is  the  discrete  Fourier  transform  of  u(x,  t). 

For  each  wavenumber,  the  Fourier  velocity  u  evolves  by 


du 

dt 


=  a[a(t)], 


(l) 
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where  a  (obtained  from  the  Navier-Stokes  equations)  represents  the  velocity  time 
derivative  in  wavenumber  space  and  is  dependent  on  u  at  all  k.  This  (vector) 
equation  is  integrated  in  time  by  an  explicit  second-order  Runga-Kutta  method. 
The  essence  of  the  pseudo -spectral  method  is  that  the  velocity  products  (that 
arise,  for  example,  in  the  convective  term)  are  evaluated  in  physical  space,  and 
then  transformed  to  wavenumber  space.  This  avoids  the  costly  evaluation  of  the 
products  in  physical  space  as  convolutions  in  wavenumber  space.  On  the  other 
hand,  spatial  derivatives  (that  would  have  to  be  approximated  in  physical  space) 
are  evaluated  without  approximation  in  wavenumber  space.  Most  of  the  compu¬ 
tational  work  (70%  of  it)  is  consumed  in  the  Fast  Fourier  Transforms  (FFT's) 
used  to  transform  the  data  between  physical  and  wavenumber  spaces. 

A  forcing  scheme  [11]  is  used  to  maintain  the  turbulence  energy  against  viscous 
decay.  This  is  achieved  by  adding  an  artificial  random  term  aF(k,  t)  to  the  right- 
hand  side  of  Eq.  (1).  It  has  been  verified  [15]  that  forcing  the  large-scale  motions 
has  a  negligible  effect  on  the  small-scale  motions  that  strain  and  bend  material 
surfaces. 

16.3.2  Implementation  on  the  IBM  3090-600E 


Rogallo's  original  pseudo-spectral  code  [14]  was  written  in  the  Vectoral  program¬ 
ming  language  and  run  on  the  Cray  1  computer  at  NASA  Ames.  Because  of  the 
very  limited  memory  of  this  computer,  a  great  deal  of  the  coding  is  concerned 
with  transferring  data  between  memory  and  secondary  storage. 

While  maintaining  the  same  numerical  algorithm,  we  completely  rewrote  the 
code  in  FORTRAN  to  exploit  fully  the  vector,  parallel  and  large-memory  capa¬ 
bilities  of  the  IBM  3090-600E.  In  a  1283  simulation  (i.e.  N  =  128)  the  grid  con¬ 
sists  of  over  2  million  nodes.  As  shown  in  Figure  149,  these  nodes  can  be 
considered  to  be  N  x-y  planes  (of  N2  nodes  each),  or  as  N  x-z  planes.  The 
essence  of  the  parallel  implementation  is  to  perform  operations  on  planes  of  data, 
assigning  each  plane  to  a  separate  process.  Thus  in  each  phase  of  the  calculations 
N  processes  are  generated,  each  pertaining  to  an  x-y  or  an  x-z  plane  of  data;  and, 
at  a  given  time,  each  of  the  6  processors  of  the  IBM  3090-600E  is  executing  one 
of  these  processes.  This  is  an  extremely  simple  and  effective  use  of  parallelism 
which  depends  upon  two  properties:  the  large  shared  memory  of  the  IBM 
3090-600E;  and,  the  ability  to  split  the  algorithm  into  operations  requiring  only 
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Figure  149.  Data  Structure.  Sketch  showing  an  N 3  of  data  as:  a)  N  x-y  planes,  and  b)  N 
x-z  planes. 


lines  or  planes  of  data.  The  speed-up  achieved  (compared  to  using  a  single 
processor)  is  greater  than  4. 

The  usual  measures  have  been  taken  to  exploit  the  speed-up  made  possible  by 
the  vector  processors  of  the  IBM  3090.  Since  about  70%  of  the  CPU  tune  is 
consumed  in  performing  FFT's,  particular  attention  was  paid  to  their  implemen¬ 
tation.  The  ESSL  vector  routines  are  used  for  this  purpose,  with  the  optimum 
data  structure  and  calling  sequence  having  been  determined  by  experimentation. 

The  computational  requirements  of  a  (128)3  simulation  are  very  large.  Each 
time  step  takes  a  total  of  2  CPU  minutes,  and  a  full  run  requires,  typically,  2,400 
time  steps.  Thus  the  total  CPU  time  required  is  about  80  hours.  However, 
because  of  the  use  of  parallelism,  the  turnaround  time  can  be  less  than  20  hours. 

16.3.3  Physical  and  Numerical  Parameters 


Results  are  reported  from  four  simulations,  which  are  distinguished  by  the 
Reynolds  number  R;.  Table  26  on  page  495  contains  a  summary  of  the  specified 
physical  and  numerical  parameters  and  the  resulting  basic  turbulence  properties. 

The  accuracy  of  the  simulations  is  determined  by  the  spatial  and  temporal  resol¬ 
ution.  Detailed  tests  have  been  performed  [11,  12,  15]  to  establish  that  resol¬ 
ution  is  excellent  for  all  the  simulations. 
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Figure  150.  Spectrum  Against  Kolmogorov-Scaled  Wavenumber.  Spectrum  k3E(k) 
(arbitrary  scale)  against  Kolmogorov-scale  wavenumber.  A  1283R^  63;  □  643R/i 
59;  dashed  line,  experimental  data  of  Comte- Bellot  &  Corrsin  [16],  65. 

To  illustrate  the  veracity  of  the  simulations,  on  Figure  150  we  show  computed 
and  experimentally  measured  dissipation  spectra.  With  E(k)  being  the  energy 
spectrum  function,  the  quantity  shown  on  the  figure  is  proportional  to  k3E(k)  — 
thus  accentuating  the  high  wavenumber  (i.e.  small  scale)  components  of  the 
velocity  field.  From  this  figure  (and  other  tests  performed  [11,  12,  15])  we  can 
conclude: 

1.  the  spatial  resolution  is  excellent  (the  whole  dissipation  spectrum  is  captured), 

2.  the  small  scale  motions  are  unaffected  by  the  details  of  the  forcing, 

3.  the  simulated  small  scales  have  the  same  statistics  as  experimentally  realized 
turbulence. 


16.4  Evolution  of  Material  Surfaces 
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16.4.1  Direct  Method 


Given  that  the  time-dependent  three-dimensional  velocity  field  u(x,  t)  is  known 
from  the  simulations  described  in  the  previous  section,  a  direct  method  of  the 
computing  the  evolution  of  a  material  surface  immediately  suggests  itself. 

For  definiteness  consider  the  material  surface  that  initially  (t  =  0)  is  the  plane 
Xi  =  0.  Let  X(t)  denote  the  location  of  a  point  on  the  surface.  Then,  since  by 
definition  the  surface  moves  with  the  fluid,  we  have: 

X(t)  =  -£-X(t)  =  u(X(t),t).  (2) 

The  direct  method  is  to  represent  the  surface  discretely  by  M  nodes,  the  m-th 
having  position  X(m)(t).  Because  turbulence  stretches  and  bends  the  surface, 
adaptive  gridding  (to  re-position  nodes  and  to  introduce  new  ones)  would  be 
required  to  maintain  resolution  of  the  surface. 

Our  results  show  that  the  deformation  of  the  surface  is  so  severe  as  to  make  this 
approach  impracticable.  During  the  simulation  the  area  of  the  surface  increases 
by  a  factor  of  1017,  and  radii  of  curvature  as  small  as  10-8  of  a  grid  spacing  occur. 
Clearly  it  is  hopeless  (and  misguided)  to  attempt  to  represent  such  a  massive 
surface  with  the  necessarily  fine  resolution. 

16.4.2  Infinitesimal  Surface  Elements 


Rather  than  attempting  to  represent  and  resolve  the  whole  surface,  we  study, 
instead,  a  large  number  (M  =  4,096  or  8,192)  of  representative  infinitesimal 
surface  elements.  A  comprehensive  account  of  the  definition  and  properties  of 
such  elements  is  provided  by  Pope  [2]. 

Figure  151  is  a  sketch  of  an  infinitesimal  surface  element  at  time  t.  Its  position 
is  X(t),  its  infinitesimal  area  is  dA(t),  and  the  unit  normal  to  the  surface  is  N(t). 
For  each  element,  a  time-dependent  Cartesian  coordinate  system  is  introduced, 
with  its  origin  at  X(t),  and  with  orthonormal  basis  vectors  e4(t)  (i  =  1,2,3).  The 
unit  vector  e3(t)  is  coincident  with  the  normal,  i.e.  e3(t)  =  N(t),  and  consequently 
ei(t)  and  e2(t)  are  in  the  tangent  plane  of  the  surface  at  X(t).  Initially  at  (t  =  0) 
ei(0)  is  specified  arbitrarily  in  the  tangent  plane,  and  £2(0)  is  determined  by 
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N  =  e 


Figure  151.  Sketch  of  an  Infinitesimal  Area  Element 


orthogonality.  Subsequently,  ei(t)  and  ^(t)  rotate  with  the  fluid.  These  specifica¬ 
tions  lead  [2]  to  the  following  evolution  equations: 


•  _j_  f  d^i  au«  i° ,  r  f  gu3  r 

**  2  dya  dyp  J  -3{  dya  j 


5U3 

dy* 


(3) 


Here  yi  and  y2  are  coordinates  in  the  ei  and  §2  directions;  Greek  suffices  take  the 
values  1  or  2;  the  summation  convention  applies;  Uj  =  ^  •  u;  and,  the  superscript 
0  indicates  that  the  quantities  are  evaluated  at  the  origin. 


With  dA(t)  being  the  infinitesimal  area  of  the  surface  element,  the  area  amplifi¬ 
cation  factor  A(t)  is  defined  by  A(t)=dA(t)/dA(0).  Note  that  A(0)  is  unity.  The 
area  increases  due  to  straining  according  to 


A  =  Aa . 


(4) 


where  a(t)= 


is  the  rate  of  strain  in  the  tangent  plane. 


The  curvature  of  the  infinitesimal  surface  element  is  completely  described  by  the 
(second-order  symmetric)  curvature  tensor  h^(t)  (see  Pope  [2]).  The  principal 
curvatures  ki(t)  and  k2(t)  are  the  eigenvalues  of  this  tensor,  with  the  convention 
ki  >  k2.  The  exact  evolution  equation  for  h^  is  [2]: 


d2V3 

dVocdyp 


(Sy^ay  Sya^/?y)’ 


(5) 
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where  %/l  2  )  dyfi  dy„  j  is  the  strain-rate  tensor  in  the  tangent  plane. 

(Note  that  a  =  sa«-) 

A  complete  set  of  equations  has  now  been  presented  for  the  surface  properties 
X,N  =  g3,  e„,  A  and  h“;5  (Eqs.  2-5).  From  given  initial  conditions,  these  can  be 
integrated,  given  the  time  series  of  u,  du.jdx,  and  cT-uJdx^  following  the  fluid 
particle.  Hence  the  curvatures  ki  and  k2  are  determined  as  the  eigenvalues  of  h^. 

Since  the  turbulence  is  isotropic,  the  initial  condition  N(0)  is  arbitrary.  For 
convenience  we  specify  gj(0)  to  be  coincident  with  the  axes  used  in  the  Direct 
Numerical  Simulation.  The  infinitesimal  surface  elements  are  specified  to  be 
plane  initially:  that  is  h^(0)  =  0. 

16.4.3  Numerical  Implementation 


As  the  turbulence  simulation  steps  through  time,  the  set  of  ordinary  differential 
equations  describing  the  infinitesimal  surface  elements  is  solved  by  a  second-order 
Runga-Kutta  method  for  each  of  the  M  elements  considered. 

The  set  of  differential  equations  for  the  elements  (Eqs.  2-5)  contains  the  velocity 
and  its  spatial  gradients  evaluated  at  the  element  location  X(t).  These  quantities 
are  obtained  from  the  velocities  at  the  grid  nodes  by  cubic-spline  interpolation, 
which  has  been  shown  to  be  extremely  accurate  [17],  The  formation  of  the 
splines  requires  significant  computational  effort.  It  is  implemented  on  the 
IBM3090  exploiting  the  parallel  processing  capabilities  in  much  the  same  way  as 
in  the  pseudo-spectral  code. 

16.4.4  Numerical  Accuracy 


The  numerical  accuracy  depends  on  the  spatial  resolution  (indicated  by  A x/»/ 
being  small),  and  on  the  temporal  resolution  (indicated  by  A f/r,  and  the  Courant 
number  C  being  small),  where  rj  and  t,  are  the  Kolmogorov  length  and  time 
scales.  For  the  Eulerian  simulation  it  has  been  demonstrated  in  previous  studies 
[11,  15,  17]  that  good  resolution  is  achieved  with  kmaxi/sTA  (corresponding  to 
Ax/>/«2)  and  C«0.5 
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The  results  show  that  surface  radii  of  curvature,  R,  less  than  a  millionth  of  the 
Kolmogorov  scale  r\  are  observed.  Since  both  Ax/R  and  the  Courant  number 
based  on  R  can  be  very  large  (  ~  105),  careful  consideration  needs  to  be  given  to 
whether  numerical  accuracy  for  the  radii  of  curvature  can,  nevertheless,  be 
claimed. 

The  sources  of  numerical  error  in  determining  ki(t)  and  k2(t)  —  beyond  those 
incurred  in  the  Eulerian  simulation  —  are  threefold.  First  there  is  the  time¬ 
stepping  error  in  integrating  the  surface -element  ordinary  differential  equations. 
Inspection  of  Eq.  (5)  suggests  that  the  time  scale  of  change  of  ka  is  no  smaller 
than  the  Kolmogorov  time  scale  rn:  this  is  confirmed  by  the  result.  Hence  the 
time  step  of  size  0.1  r7  is  sufficiently  small,  as  tests  verify. 

Second,  there  is  some  error  involved  in  interpolating  for  the  velocity  derivatives. 
As  the  tests  performed  by  Yeung  and  Pope  [17]  show,  with  the  current  resol¬ 
ution  (kmax*/«1.5)  and  using  cubic  spline  interpolation,  this  error  is  less  than  1%. 

Third,  there  is  a  numerical  error  (again  investigated  by  Yeung  and  Pope  [17])  in 
the  integration  of  Eq.  (2)  to  determine  the  element  location  X(t).  Given  the 
dispersive  nature  of  turbulence,  for  t  large  compared  to  Tn,  this  error  could  be 
large  —  certainly  large  compared  to  10-V  But  the  error  in  X(t)  can,  alterna¬ 
tively,  be  viewed  as  a  small  error  in  the  initial  condition.  That  is,  the  numerically 
determined  particle  position  X(t)  is  the  exact  position  of  the  particle  originating 
from  X(0)  +  <5X  where  |<5X|  is  small  (compared  to  rj  ).  And  the  time  series  of 
the  velocity  gradients  of  the  particles  originating  from  X(0)  and  X(0)  +  <5X  differ 
little.  Since  we  are  interested  in  the  statistics  of  a  statistically  homogeneous 
surface,  the  precise  initial  condition  of  the  surface  elements  considered  is  unim¬ 
portant. 

In  summary,  in  the  current  method  of  tracking  infinitesimal  surface  elements, 
good  resolution  of  the  Kolmogorov  length  and  time  scales  is  sufficient  for  the 
accurate  calculation  of  surface  statistics.  Resolution  on  the  scales  of  the  surface 
radii  or  curvature  is  not  required.  This  conclusion  is  in  marked  contrast  to  that 
for  the  direct  method  (16.4.1,  “Direct  Method”  on  page  480)  in  which  the  whole 
surface  is  represented  numerically. 
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Figure  152.  Probability  Density  Functions  (pdf)  of  Surface  Strain.  Probability 
density  functions  of  Kolmogorov-scaled  strains  acting  on  material  and  fixed  surfaces. 
Ri  93  :  A  a* (area- weighted);  □  a*  (unweighted);  O  a(for  fixed  surface). 
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16.5  Results  and  Discussion 

16.5.1  Straining  on  Material  Surfaces 


In  this  subsection  we  first  present  statistics  characterizing  the  straining  of  material 
surfaces.  The  fundamental  quantity  is  a(t),  the  (total)  rate  of  strain  in  the  tangent 
plane  of  the  material  surface,  which  determines  the  rate  of  surface  area  increase 
(Eq.  4).  It  is  appropriate  to  consider  one-time  area-weighted  statistics,  such  that 
equal  areas  contribute  equally  to  the  mean.  We  note  that  straining  statistics  from 
our  four  simulations  display  no  significant  Reynolds  number  dependence  [3]. 

The  rate  of  strain  in  the  tangent  plane  of  a  slowly  propagating  surface  (w  <  v,, 
where  v,  is  the  Kolmogorov  velocity  scale)  is  well  approximated  by  a;  whereas  for 
a  rapidly  propagating  surface  (w>v^ )  it  is  close  to  a,  the  rate  of  strain  acting  on  a 
randomly  oriented  element,  which  is  statistically  identical  to  —  dui/dxi. 
Figure  152  shows  the  area- weighted  (f£(a))  and  unweighted  (fQ(tf))  probability 
density  functions  (pdfs)  of  a,  and  the  pdf  of  a(f.(a)«f2f(a)).  The  strain  rates  are 
Kolmogorov-scaled  (a*  =  aT,,;a*  =  ar,). 

It  is  seen  that,  relative  to  fa(a),  f£(a)  is  displaced  to  the  right,  since  large  areas  are 
associated  with  strong  stretching  (positive  a).  In  fact,  the  area-weighted  mean 
<  a*  >A  =  0.28  is  almost  twice  the  unweighted  mean.  The  mean  of  a  is  iden¬ 
tically  zero,  and  its  skewness  is  nearly  0.5,  consistent  with  previous  results  [18]. 
The  variance  of  a*  is  close  to  1/15,  its  theoretical  value  in  isotropic  turbulence. 

Two  principal  strain  rates,  Si  and  S2  (with  Si  >  S2,  and  Si  +  S2  =  a),  are  the 
eigenvalues  of  the  strain-rate  tensor  in  the  tangent  plane  of  a  material  surface  sap. 
The  area-weighted  joint  pdf  of  Kolmogorov- scaled  Si  and  S2  (Sf  =  Sit,,  Sf  =  S2t„) 
is  shown  in  Figure  153.  The  isoprobability  contours  indicate  that  Si  is  almost 
always  positive  (with  98%  probability),  and  that  there  is  relatively  high  proba¬ 
bility  for  the  magnitude  of  S2  to  be  small.  The  mean  values  of  Sf  and  Sf  are 
found  to  be  approximately  0.315  and  -0.035  respectively.  The  probability  of  net 
compressive  straining  in  the  tangent  plane,  i.e.,  Si  +  S2  <  0,  is  about  19%. 

We  now  examine  the  conjectures  made  by  Batchelor  [1,  7]  on  the  straining  of 
material  surfaces  (see  16.2,  “Introduction”).  As  noted  above,  the  area-weighted 
mean  value  of  the  surface  strain  is  approximately  <  a  >A  =  0.28/r,.  Thus  the 
time  scale  of  surface  strain  can  be  taken  as  l/<  a  >a^3.5t7.  The  question  of 
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Figure  153.  Joint  pdf  of  Principal  Surface  Strains.  Contour  plot  of  the  area-weighted 
joint  pdf  of  the  Kolmogorov-scaled  principal  surface  strain  rates  S*i  and  S*2.  The 
indicated  regions  and  their  probabilities  are: 

A:  Si  >  S2  >  0  >  SN,  Prob  (A)  =  0.498 
B:  Si  >  0  >  S2  >  SN,  Prob  (B)  =  0.201 
C:  Si  >  0  >  SN  >  S 2 ,  Prob  (C)  =  0.110 
D:  Si  >  SN  >  0  >  S2,  Prob  (D)  =  0.100 
E:  SN  >  Si  >  0  >  S2>  Prob  (E)  =  0.077 
F:  SN  >  0  >  Si  >  S2,  Prob  (F)  =  0.014 
where  S^  =  -a  =  -  Si  -  S2 

whether  the  straining  is  persistent  is  best  addressed  by  examining  the  two-time 
statistics  of  a,  in  particular  its  autocorrelation  function.  This  is  shown  in 
Figure  154  for  the  38  and  93  cases.  It  may  be  seen  that  the  autocorrelation 
decays  to  small  values  well  before  a  time  lag  of  3.5r^  is  reached,  and  in  fact 
crosses  the  zero  line  at  about  2.75rr  The  integral  time  scale  of  a  is  about 
1.0t^(0.96t^  at  R^  38,  1.03t^  at  R;  93),  which  is  small  compared  to  3.5t7.  From 
these  observations,  we  conclude  that  the  time  scale  of  change  of  strain  is  certainly 
not  large  compared  to  the  time  scale  of  strain  itself  —  i.e.,  the  straining  on  a 
material  surface  is  not  persistent. 
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Figure  154.  Autocorrelations  of  Surface  Strain  Rate  a(t).  Vertical  axis  is  /?«(?*).  Hor¬ 
izontal  axis  is  t*  =  t/t7.  A  =  38  O  =  93. 

The  second  conjecture  mentioned  in  16.2,  “Introduction”  is  that  the  surface 
element  becomes  aligned  with  the  principal  axes  corresponding  to  the  two 
greatest  principal  strain  rates.  Let  TA  be  the  angle  between  the  normal  to  the 
surface  N,  and  the  principal  axis  corresponding  to  the  least  principal  strain  rate. 
Then  the  conjecture  is  equivalent  to  TA  =  0.  Figure  155  shows  the  evolution  of 
the  mean  <TA>  from  the  simulations.  Initially  <TA>  has  a  value  close  to 
unity,  corresponding  to  the  initial  condition  of  TA  being  uniformly  random. 
After  a  transient,  <  TA  >  attains  a  value  close  to  0.75  radians.  Clearly  then,  the 
conjecture  (FA#0)  is  far  from  the  truth. 
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Figure  155.  Evolution  of  The  Mean  Angle  <rA>.  The  evolution  of  the  mean  angle  < 
Ta>  between  the  normal  to  the  surface  and  the  principal  axis  corresponding  to  the 
least  principal  strain  rate:  (O),  38;  (•),  63;  (★),  R^  90. 

16.5.2  Curvature  of  Material  Surfaces 


Curvature  statistics  have  been  determined  for  the  38  simulation,  by  calculating 
the  properties  of  an  ensemble  of  8,192  infinitesimal  surface  elements.  For  each 
element  the  principal  curvatures  ki  and  k2  are  determined  as  the  eigenvalues  of 
h|^;  and  the  area  amplification  A  is  used  to  construct  area-weighted  statistics  [3], 

The  curvature  is  characterized  by  the  following  quantities:  the  normalized 
mean- square  curvature,  M*(f),  the  logarithm  of  curvature,  L,  and  ,  the  normal¬ 
ized  radius  of  curvature  of  the  material  element  R*(7): 

L{t)  =  In  M*(t)  =  y  \t)  +  kl(t)),  R *  = 


(6) 
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Figure  156.  Standardized  pdf  of  Logarithm  of  Curvature.  Vertical  axis  is  logi0[i^(/ )]. 

A 

Horizontal  axis  is  / .  Standardized  area- weighted  pdf  of  the  logarithm  of  curvature 
L  =  InM*  (Eq.  12).  Dashed  line  corresponds  to  lognormal  distribution  for  M*; 
straight  line  corresponds  to  Eq.  (16). 

It  is  found  that  an  essentially  statistically  stationary  state  is  attained  after  about 
15  Kolmogorov  time  scales.  (Some  subtleties  about  this  state,  and  methods  of 
reducing  statistical  errors  in  the  state  are  discussed  in  Yeung  et  al.[3]  and  Pope  et 
al.[12]).  The  area- weighted  pdf  pA(i)  of  the  standardized  logarithm  of  curvature, 
L  =  (L  —  <L  >a)/<7a,  is  shown  in  Figure  156.  The  mean  and  standard  deviation 
of  L  are  <  L  >A  =  -4.77  and  oA  =  2.18.  The  dashed  line  on  the  figure  is  the 
parabola  corresponding  to  the  pdf  if  M*  were  log-normally  distributed.  Clearly 
the  log-normal  distribution  does  not,  even  qualitatively,  describe  the  shape  of  the 
pdf. 

It  appears  from  the  figure,  that  for  large  curvatures,  the  pdf  of  L,  pA(i),  has  the 
asymptotic  form 

p  (i)~be 


(7) 
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Figure  157.  Time  Series  of  Normalized  rms  Curvature.  Verticle  axis  is  ^ jM\t )  .  Hori¬ 
zontal  axis  is  r/T^.  Normalized  rms  curvature  against  time  for  a  typical  element  (O), 
and  for  the  element  that  attains  the  greatest  curvature  ( A ). 

with  b  =  0.0172  and  c  =  0.55.  An  immediate  consequence  of  this  result  is  that 
the  expectation  of  the  mean-square  curvature  <  M*  >A  is  infinite.  For,  if  the 
integral  converged,  <  M*  >A  would  be  given  by 

pOO 

<  M*  >A  =  I  e  pA(i)di.  (8) 

— GO 

Clearly,  with  the  asymptotic  form  of  Eq.  (7)  with  c  <  1,  the  integral  does  not 
converge. 

To  illustrate  directly  the  occurrence  of  large  curvatures,  in  Figure  157  we  show 
the  time  series  Jm*  for  the  surface  element  that  attains  the  largest  value  of  M* 
It  can  be  seen  that  at  parts  it  experiences  a  rapid  rise  in  M* ,  at  the  approximate 
exponential  rate  of  dL/dt»1.8/r,.  The  peak  value  of  M*  corresponds  to  a  radius 
of  curvature  R  less  than  10~8^. 
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We  examine  in  Figure  1 58  the  area- weighted  pdf  pR(r)  and  the  unweighted  pdf 
pR(r)  of  the  normalized  mean  radius  of  curvature  (Eq.  6).  The  mean  <  R*  >A  is 
12.0  (corresponding  to  about  half  the  integral  length  scale).  Even  though  the 
mean-square  curvature  <  M*  >A  tends  to  infinity,  Figure  1 58  shows  that  only 
about  5%  of  the  surface  has  mean  radius  of  curvature  R  smaller  than  the 
Kolmogorov  length  scale  r\.  The  unweighted  and  the  area- weighted  pdfs  differ 
greatly,  especially  at  the  origin.  The  implication  is  that  highly  curved  elements 
tend  to  have  less  area  than  more  mildly  curved  elements. 

The  shape  of  the  surface  element  at  each  point  is  determined  by  the  relative 
values  of  the  two  principal  curvatures.  We  define  the  shape  parameter  9  by 
9  =  ks/k„  where,  ks  and  kt  are  the  smaller  and  larger  of  /^andfe  in  absolute  magni¬ 
tude.  Possible  values  of  6  lie  between  —1  and  4-1.  The  value  9  —  1  corresponds 
to  a  spherical  element;  the  value  6  =  0  corresponds  to  a  cylindrical  element;  and 
the  value  0  =  —  1  corresponds  to  a  pseudo -spherical  element. 


Figure  158.  PDF  of  Normalized  Mean  Radius  of  Curvature.  Area-weighted  (O)  and 
unweighted  ( A )  pdf  of  normalized  mean  radius  of  curvature. 
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Figure  159.  Area  Pdf  of  the  Shape  Parameter  9.  Area  weighted  pdf  of  the  shape 
parameter  ©. 

In  Figure  159  we  present  the  pdf  of  0 ,  and  in  Figure  160  a  contour-plot  of  the 
joint  pdf  of  6  and  L  is  presented.  From  Figure  159  it  is  clear  that  the  cylindrical 
shape  is  more  probable,  and  the  probability  of  a  material  element  being  spherical 
in  shape  is  very  small  compared  to  either  the  cylindrical  shape  or  the  pseudo- 
spherical  shape.  The  joint  pdf  of  9  and  L  in  Figure  160  shows  that  cylindrical 
shaped  elements  are  associated  with  much  higher  than  average  curvatures, 
whereas  the  pseudo-spherical  elements  are  more  mildly  curved. 


16.6  Conclusions 


Since  the  earliest  work  on  turbulence  almost  a  century  ago,  there  has  been  no 
shortage  of  imaginative  theories.  But  these  theories  have  been  based  on  hypoth¬ 
eses  and  conjectures  that,  at  the  time,  were  impossible  to  test.  The  present  work 
clearly  demonstrates  that  turbulence  research  is  entering  a  new  era,  one  in  which 
supercomputers  are  used  to  extract  detailed  information  about  turbulence  — 
information  that  can  be  used  to  test  hypotheses  and  suggest  new  ones. 
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Contrary  to  previous  wisdom,  it  has  been  shown  that  turbulent  straining  is 
fleeting,  rather  than  persistent.  As  a  consequence  material  elements  do  not 
become  aligned  with  the  principal  axes  of  the  strain  time. 

It  is  found  that  material  surfaces  become  extremely  contorted  due  to  the 
straining  and  bending  of  the  turbulence.  Radii  of  curvature  as  small  as  10~8  of 


e 


Figure  160.  Joint  pdf  of  In  M*  and  0.  Contour  plot  of  the  (unweighted)  joint  pdf  of  L  - 
InM*  and  0. 
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the  Kolmogorov  scale  are  observed.  And  these  highly-curved  elements  are  almost 
cylindrical  in  shape. 

The  simulations  described  here  are  for  the  simplest  possible  turbulent  flow  — 
homogeneous  isotropic  turbulence.  But  the  result  have  general  applicability  since 
the  behavior  of  the  small  scales  is  (nearly)  universal  —  as  evidenced  by  the  lack  of 
dependence  of  our  results  of  Reynolds  number  R^. 
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Table  26.  Numerical  Parameters  and  Eulerian  Statistics 

Taylor- scale 

Reynolds 

number 

Ra 

38 

59 

63 

90 

93 

Grid  Size 

N 

64 

64 

128 

128 

128 

Length  of  sol¬ 
ution  domain 

Lo 

2  7T 

2  7t 

2  7T 

2  n 

2  n 

Kinematic 

viscosity 

V 

0.025 

0.0158 

0.0105 

0.006546 

0.006546 

Turbulence 

intensity 

u' 

1.60 

1.276 

1.637 

1.274 

1.356 

Dissipation  rate 

<£> 

2.69 

0.775 

2.673 

0.780 

0.893 

Longitudinal 
integral  length 
scale  L\ 

Li 

(T^> 

0.369 

0.528 

0.321 

0.448 

0.398 

Dissipation 
time  scale 

t£  =  —  w'2/  <  £  > 

1.43 

3.184 

1.510 

3.174 

3.099 

Eddy  turnover 
time  Te  =  Li/u' 

T./t, 

0.507 

0.406 

0.407 

0.343 

0.297 

Kolmogorov 
time  scale  'tn 

0.067 

0.045 

0.041 

0.029 

0.0028 

Duration  of 
simulation  T 

T/Te 

13.8 

8.35 

5.85 

5.51 

6.53 

Time  step  At 

Af/T, 

0.52 

0.042 

0.024 

0.027 

0.029 

Kolmogorov 
length  scale  r\ 

n!L  i 

0.042 

0.029 

0.025 

0.018 

0.019 

Maximum 

resolved 

wavenumber 

kmax 

kmax^? 

1.48 

1.43 

1.54 

1.48 

1.42 

Taylor  micro¬ 
scale  A 

ML  1 

0.521 

0.437 

0.399 

0.326 

0.359 

496 


Pope,  Yeung  and  Girimaji 


16.8  References 

1.  G.K.  Batchelor,  J.  Fluid  Mech.,  5,  113  (1959). 

2.  S.B.  Pope,  Int.  J.  Engr'ng  Sci.,  26,  445  (1988). 

3.  P.K.  Yeung,  S.S.  Girimaji  &  S.B.  Pope,  Combust.  Flame ,  79,  340  (1990). 

4.  S.B.  Pope,  Ann.  Rev.  Fluid  Mech ,  19,  237  (1987). 

5.  B.  Karlovitz,  D.W.  Denniston,  Jr.,  D.M.  Knappschaefer  and  F.E.  Wells,  Fourth  Symp.  (lnt’l) 
on  Combust.,  p.  613,  Williams  &  Wilkins  (1953). 

6.  M.  Matalon  &  B.J.  Matkowsky,  J.  Fluid  Mech.,  124,  239  (1982). 

7.  G.K.  Batchelor,  Proc.  Roy.  Soc.,  A,  213,  349  (1952). 

8.  A.S.  Monin  &  A.M.  Yaglom,  Statistical  Fluid  Mechanics,  Vol.  2  (J.L.  Lumley,  Ed.)  M.I.T  Press 
(1975). 

9.  S.  Corrsin,  J.  Atmos.  Sci.,  20,  115  (1962). 

10.  A.M.  Klimov,  Sov.  Phys.  Doki,  20,  168  (1975). 

11.  P.K.  Yeung  &  S.B.  Pope,  J.  Fluid  Mech.,  207,  531  (1989). 

12.  S.B.  Pope,  P.K.  Yeung  and  S.S.  Girimaji,  Phys.  Fluids,  A,  1  2010  (1989). 

13.  S.S.  Girimaji  &  S.B.  Pope,  “Material  Element  Deformation  in  Isotropic  Turbulence,”  J.  Fluid 
Mech.  (1990)(in  press). 

14.  R.S.  Rogallo,  “Numerical  Experiments  in  Homogeneous  Turbulence,”  NASA  TM  81315  (1981). 

15.  V.  Eswaran  &  S.B.  Pope,  Comput.  Fluids,  16,  157  (1988). 

16.  G.  Comte- Bellot  and  S.  Corrsin,  J.  Fluid  Mech,  48,  273  (1971). 

17.  P.K.  Yeung  &  S.B.  Pope,  J.  Comput.  Phys.,  19,  373  (1988). 

18.  R.M.  Kerr,  J.  Fluid  Mech,  153,  31  (1985). 


